% DO THIS BEFORE EXECUTING
% addpath(genpath('/media/sda2/Users/Nermine/Desktop/code/mpt'));
%addpath(genpath('/media/sda2/Users/Nermine/Desktop/code/sdp/gurobi301/gurobi_mex_v1.20'));
% Syntax:
%     x = gurobi_mex(c, objtype, A, b, contypes, lb, ub, vartypes); 
%     x = gurobi_mex(c, objtype, A, b, contypes, lb, ub, vartypes, options); 
%     [x,val] = gurobi_mex(...); 
%     [x,val,flag] = gurobi_mex(...); 
%     [x,val,flag,output] = gurobi_mex(...); 
%     [x,val,flag,output,lambda] = gurobi_mex(...); 
%     *  c: objective coefficient vector, double. 
% 
%     [] (empty array) means uniformly 0 coefficients, and scalar means all coefficients equal to scalar. 
% 
%     * objtype: 1 (minimization) or -1 (maximization).
%     * A: constraint coefficient matrix, double, sparse.
%     * b: constraint right-hand side vector, double. 
% 
%     Gurobi takes a dense vector for this input. If a sparse vector is specified, it is converted to full by Gurobi Mex. 
% 
%     * contypes: constraint types. Char array of '>', '<', '='. 
% 
%     Warning: '>=' means two constraints instead of one an inequality constraint. 
%     Example: '>><=' means the first two constraints have greater or equal to signs, the third has less than or equal to sign, and the last is an equality constraint. 
%     If a single character is specified, all constraints are uniformly signed to the corresponding type. 
% 
%     * lb: variable lower bound vector, double. 
% 
%     [] (empty array) means 0 lower bound. -inf means no lower bound. scalar means a uniform lower bound equal to scalar. 
% 
%     * ub: variable upper bound vector, double. 
% 
%     [] (empty array) means no (or infinity) upper bound. scalar means a uniform upper bound equal to scalar. 
% 
%     * vartypes: variable types. Char array of chars 'C', 'B', 'I', 'S', 'N'. C for continuous; B for binary; I for integer; S for semi-continuous; N for semi-integer. [] (empty array) means all variables are continuous. 
% 
%     Example: 'CCCCC' stands for five continuous variables. 
%     Note that semi-continuous variables are variables that must take a value between their minimum and maximum or zero. Semi-integer variables are similarly defined. 
%     If a single character is specified, all variables are uniformly signed to the corresponding type. 
% 
%     * options: optional structure that may contain one or more of the following fields: (see Gurobi's parameter help for their allowed values. Also, see examples below.)
%           o options.IterationLimit: see Gurobi's parameter help.
%           o options.FeasibilityTol: see Gurobi's parameter help.
%           o options.IntFeasTol: see Gurobi's parameter help.
%           o options.OptimalityTol: see Gurobi's parameter help.
%           o options.MIPGap: see Gurobi's parameter help.
%           o options.LPMethod: see Gurobi's parameter help.
%           o options.Presolve: see Gurobi's parameter help.
%           o options.TimeLimit: see Gurobi's parameter help.
%           o options.Threads: see Gurobi's parameter help.
%           o options.DisplayInterval: Gurobi's Callback screen output interval. See Gurobi's parameter help.   0 means no Gurobi message.
%           o options.Display: Gurobi Mex's screen output level. 0 for no output; 1 for error only; 2 (default) for normal output. For complete silence, both DisplayInterval and Display need to be set to 0.
%           o options.LogFile: char array of the name of log file. options.UseLogfile is no longer used.
%           o options.WriteToFile: char array of the name of the file to which optimization data is written. See Gurobi C-Reference entry GRBwrite for supported formats. This option helps one verify whether the model is correctly passed to Gurobi. 
% 
% Output Description
% 
%     * x: primal solution vector; empty if Gurobi encounters errors or stops early (in this case, check output flag).
%     * val: optimal objective value; empty if Gurobi encounters errors or stops early.
%     * flag: value meanings:
%           o 1 for not started
%           o 2 for optimal
%           o 3 for infeasible
%           o 4 for infeasible or unbounded
%           o 5 for unbounded
%           o 6 for objective worse than user-specified cutoff
%           o 7 for reaching iteration limit
%           o 8 for reaching node limit
%           o 9 for reaching time limit
%           o 10 for reaching solution limit
%           o 11 for user interruption
%           o 12 for numerical difficulties
%           o 13 for suboptimal solution (Gurobi 3 and later) 
%     * output: structure contains the following fields
%           o output.IterCount: number of Simplex iterations
%           o output.Runtime: running time in seconds
%           o output.ErrorMsg: contains Gurobi error message, if any 
%     * lambda: Lagrange multipliers. Because solving MIPs gives no such output, do not ask for this output for MIPs. 


clc;
clear;

N = 6;
M = 100000;

%weights
q = [0.1 0.4 0.5 0.1 0.2 0.5]; % randomly generated
%q = [1 1 1 1 1 1];

C = [0 12  M  M  9 16; %Hardcoded N here
    12  0 19 12  M 15;
     M 19  0 21  M 17; 
     M 12 21  0 10 16;
     9  M  M 10  0 10;
    16 15 17 16 10  0];

c = [reshape(C,N^2,1);zeros(N^2,1)]; % vecotirzed form of CMatrix is required as input to gurobi

objtype = 1;              % 1 for minimize, -1 for maximize


% Number of variables: integers N^2 for xij and N^2 binary for yij
% Creating the R matrix required for one of the constraints

R = (sum(q)-q(1))*ones(N,N);
for i=1:N
    for j=1:N
        if(j==1)
            R(i,j) = q(1);
        end
        if(i==1)
            R(i,j) = sum(q);
        end
    end
end



A00 = [diag([1 zeros(1,N-1)]) zeros(N)];
A01 = [zeros(N) diag([1 zeros(1,N-1)])];
for i=1:N
 A0(i,:) = reshape(circshift(A00,[i-1,i-1]),2*N^2,1)'; 
 A0(i+N,:) = reshape(circshift(A01,[i-1,i-1]),2*N^2,1)';
end

A2(1,:) = [zeros(1,N^2) ones(1,N) zeros(1,N^2-N)]; %colsum of yij
A1(1,:) = [zeros(1,N^2) reshape([ones(N,1) zeros(N,N-1)]',1,N^2)]; %rowsum of yij
for i=2:N
    A2(i,:) = circshift(A2(1,:)',N*(i-1))';
    A1(i,:) = circshift(A1(1,:)',i-1)';
end

A3 = [ones(N,1); zeros(N^2-N,1); zeros(N^2,1)]'; % the Nth leg flow value is 1 back to the first node. column sum


A42(1,:) = [ones(1,N) zeros(1,N^2-N) zeros(1,N^2)]; %colsum of xij
A41(1,:) = [reshape([ones(N,1) zeros(N,N-1)]',1,N^2) zeros(1,N^2)]; %rowsum of xij
A4(1,:) = A42(1,:) - A41(1,:);
for i=2:N
    A42(i,:) = circshift(A42(1,:)',N*(i-1))';
    A41(i,:) = circshift(A41(1,:)',i-1)';
    A4(i,:) = A42(i,:)-A41(i,:);
end
bA4 = q';
bA4(1) = bA4(1) - sum(q);

A5 = zeros(N^2,2*N^2);
for i=1:N
    for j=1:N
        A50 = [zeros(N) zeros(N)];
        A50(i,j) = 1; % for xij
        A50(i,j+N) = -R(i,j); % for yij
        A5(N*(i-1)+j,:) = reshape(A50,2*N^2,1)';
    end
end

randomA6 = [reshape([ones(1,N); zeros(N-1,N)],1,N^2) zeros(1,N^2)]; % the 1st leg flow value is N exiting from first node. row sum

A =  sparse([A0; A1; A2; A3; A4; A5]);
b = [zeros(2*N,1); ones(N,1);ones(N,1); q(1); bA4; zeros(N^2,1)];
contypes = '===============================<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<';


lb = zeros(2*N^2,1); % scalar means a uniform lower bound equal to scalar (which is zero here)
ub = [sum(q)*ones(N^2,1);ones(N^2,1)]; % using loosely somewhat. Shoudl Rij figure here?
vtypes = 'CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB';

clear opts
opts.IterationLimit = 2000;
opts.FeasibilityTol = 1e-6;
opts.IntFeasTol = 1e-5;
opts.OptimalityTol = 1e-6;
opts.LPMethod = 1;         % 0 - primal, 1 - dual
opts.Presolve = -1;        % -1 - auto, 0 - no, 1 - conserv, 2 - aggressive
opts.Display = 1;
opts.LogFile = 'weighted_fischetti_gurobi_mex_MIP.log';
opts.WriteToFile = 'weighted_fischetti_gurobi_mex_MIP.mps';

[x,val,exitflag,output] = gurobi_mex(c,objtype,A,b,contypes,lb,ub,vtypes,opts);

reshape(x(1:36),6,6)